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Noise reduction for flows using nonlinear constraints 
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On the basis of a local-projective with nonhnear constraints (LPNC) approach (see K. Urbanowicz, 
J. A. Holyst, T. Stemler and H. Benner, Acta Phys. Pol B 35 (9), 2175, 2004) we develop a method of 
noise reduction in time series that makes use of constraints appearing due to the continuous character 
of flows. As opposed to local-projective methods in our method we do not need to determine the 
Jacobi matrix. The approach has been successfully applied for separating a signal from noise in 
the Lorenz model and in noisy experimental data obtained from an electronic Chua circuit. The 
method was then applied for filtering noise in human voice. 

PACS numbers: 05.45.Tp,05.40.Ca 



I. INTRODUCTION 



It is common that observed data are contaminated by 
noise (for a review of methods of nonhnear time series 
analysis see 0, 0; 0)- The presence of noise can sub- 
stantially affect such system parameters as dimension, 
entropy or Lyapunov exponents 0. In fact noise can 
completelv destroy the fractal structure of a chaotic at- 
tractor and even 2% of noise can make a dimension 
calculation misleading It follows that both from the 
theoretical as well as from the practical point of view it 
is desirable to reduce the noise level. Thanks to noise 
reduction |1 i, i, i, [13 [IE [11 [11 El [11 [11 El El [13 
it is possible e.g. to restore the hidden structure of an 
attractor which is smeared out by noise, as well as to 
improve the quality of predictions. 

Every method of noise reduction assumes that it is 
possible to distinguish between noise and a clean sig- 
nal on the basis of some objective criteria. Conventional 
methods such as linear filters use a power spectrum for 
this purpose. Low pass filters assume that a clean signal 
has some typical low frequency, respectively it is true for 
high pass filters. It follows that these methods are con- 
venient for a regular source which generates a periodic 
or a quasi-periodic signal. In the case of chaotic signals 
linear filters cannot be used for noise reduction without 
a substantial disturbance of the clean signal. The reason 
is the broad-band spectrum of chaotic signals. It follows 
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that for chaotic systems we make use of another generic 
feature of dissipative motion that is located on attractors 
consisting of subset of smooth manifolds of an admissi- 
ble phase space. As result corresponding state vectors 
reconstructed from time delay variables are limited to 
geometric objects that can be locally linearized. This 
fact is a common background of all local projective (LP) 
methods of noise reduction. 

Besides the LP approach there are also noise reduc- 
tion methods that approximate an unknown equation of 
motion and use it to find corrections to state vectors. 
Such methods make use of neural networks E2 '^^ ^ 
genetic programming Ell one has to assume some 
basis functions e.g. radial basis functions to recon- 
struct the equation of motion. Another group of methods 
are modified linear filters e.g. the Wiener filter El|) the 
Kalman filter El| > or methods based on wavelet analysis 
E3- Applications of these methods are limited to sys- 
tems with large sampling frequencies, and they are con- 
fined to the locally linear nearest neighborhood of every 
point in phase space. 

The method described in this paper can be considered 
as an extension of a local-projective with nonlinear con- 
straints (LPNC) approach that was introduced in Ref. Q]. 
We call our method the local projection with nonlinear 
constraints for flows (LPNCF). The method takes into 
account natural constraints that occur due to the contin- 
uous behavior of flows. 

The paper is organized as follows. In the following sec- 
tion we shall present the LPNCF method and the com- 
parison with LP methods is show in Sec. IIIII In Sec. IIVI 
we present examples of noise reduction and application 
to human voice filtering. 
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II. THE LPNCF METHOD 

In Ref Pi the LPNC method of noise reduction of de- 
terministic signal is presented. In this paper we intro- 
duce a method that is based on the formulation given in 
Ref 0] but it brings much better results as compared to 
LP approach. 

Let {xi} for i = 1, 2, . . . , be the time series. The 
corresponding clean signal we denote as {ii}, so when 
the measurement noise {rji} is present we come to the 
formula Xi = Xi + rji for i = 1, 2, . . . , iV. We can define 
the time delay vectors Xi = {xi^Xi-r, - ■ ■ ,Xi^[d-i)T) as 
our points in the reconstructed phase space. Then we 
can find two nearest neighbors Xfe,Xj e to vector 

x„ (X^^ is the set of nearest neighborhood of the point 
x„). Let us introduce the following function P] 

G„(s) = Xn-s {Xk+l-s - Xj+i-s) 

for s = 0, — 1. The function G„(s) vanishes for 

clean one-dimensional systems because it appears as a 
constraint after eliminating a and b from the following 
equations: 

Xn+l = aXn + b 

Xk+1 = axk + b 

Xj+i — axj -\- b. (2) 

In the case of higher dimensional systems the function 
G„(s) does not always vanish but is altering slowly in 
time for dense sampling. This is because the absolut 
value of the term G„{s) is a function of difference of 
neighboring data {xk+i~s — a^j+i-s) etc., which evolve 
smoothly in time (near neighbors behave similar in con- 
secutive time steps) . Now one can check that for a highly 
sampled clean dynamics there can be derived such a con- 
straint 

m — 1 
fc=0 

int{log2{k)) 

{l = k+ int{k/2')) (3) 

where int{z) is a integer part of z and log2{z) is a loga- 
rithm with a base 2 from z. For example with m = 8 the 



formula (O gives the following 

Cl - G„(0) - G„(l) - G„(2) + G„(3) 

-G„(4) + G„(5) + G„(6)-G„(7). (4) 

Such a criterium for a constraint can be understood easily 
if we notice that all elements G„(s) have almost the same 
value for clean data and small s. Using this we force 
the second element of constraint I^J G„(l) to take the 
oppose sign as the first element G„(0). Then the group 
consisting of the third G„(2) and fourth G„(3) elements 
should have the oppose sign to the group of the first and 
second element of the constraint Q). If we know that 
elements Gn(s) are slightly changing with increasing s 
the constraint Q should vanish for clean data and for 
large enough m. 

Similarly as in LP methods the constraints © are 
ensured in this approach by application of the method 
of Lagrange multipliers to an appropriate cost function. 
Since we expect that corrections to noisy data should be 
as small as possible, the cost function is assumed to be 
the sum of squared corrections "S* = X^fLi i^^s)'^ ■ 

It follows that we are looking for the minimum of the 
functional 

N N 
n— 1 n— 1 

After finding zero points of 2N partial derivatives one 
gets 2N equations with 2N unknown variables Sxn and 
A„. However, in such a case the derivatives of the func- 
tional (O are nonlinear functions of these variables. For 
simplicity of computing we are interested to pose our 
problem in such a way that linear equations appear which 
can be solved by standard matrix algebra. To understand 
the role of nonlinearity let us write the terms G„(s) in 
constraint C™ in such a way that an explicit dependence 
on the unknown variables is seen 

G„(s) = G (X„_s, X„_s+i) + G (<5X„_s, X„_s+i) + 

G {Xn-s,SXn-s+l) + G {SXn-s,SXn-s+l) ■ (6) 

Here we introduced the following notation 



(Xj^_s , X-^i_s^l ) — Xji^g (^Xk — s-\-l -^j — s+l^ '^'^k—s {-^j — s-l-l ^n— s+l) Xj — g (^Xji — g-^i X/^— s+l) 
G {5Xn-s,Xn-s+l) = SXn-s (xk-s+l ~ Xj-s+l) + Sxk-s {Xj-s+1 ~ Xn-s+l) + Sxj^g (Xn-s+l " Xk-s+l) 
G (X„_s, SXn-s+l) = Xn-s (Sxk-s+l " Sxj-s+l) + Xk-s {SXj^s+1 ^ SXn-s+l) + Xj-s {SXn-s+1 " Sxk-s+l) 
G {6Xn-s, SXn-s+l) = SXn-s (Sxk-s+l ^ SXj-s+l) + Sxk-s (fej-s+l ^ SXn-s+l) + SXj-s (SXn-s+l " Sxk-s+l) , (7) 
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where X„_s = {xn-s,Xk-s,Xj-s}, SXn-s = 
{Sxn-s, Sxk-s, Sxj^s}, E^nd Xfc,Xj are the near- 
est neighbors of x„. Indices are defined as 
{n, J, A: : x„, Xfc, Xj G X^^}. In the case of uncor- 
related noise and under the assumption that the 
introduced corrections completely reduce the noise effect 
dxs = —rjs i^s=i,...,N) one can neglect the nonlinear 
terms in Eqs. Q i.e. 

m 

Y,G{S^n-s,SXn-s+i)=0 yn = l,...,N. (8) 

s=0 

In the equation (|SJ) we use the fact that {iji) = and 

Taking into account the assumption one can write 
the following linear equation for the problem Q 

M ■ (5X = B, (9) 

where M is a matrix containing constant ele- 
ments, B is a constant vector, and JX^ = 
{6xi, Sx2, ... ,SxN, ^1,^2, ■■■ , Aat) are vector dependent 
variables (T - transposition). In practice it is very diffi- 
cult or even impossible to find the solution of the equa- 
tion ((nj for large N. First, it is time consuming to solve 
a linear equation with a matrix 2N x 2N matrix for 
A'' > 1000. Second, when M becomes singular the es- 
timation error of the inverse matrix is very large. 
Third, we cannot always find the true nearest neighbors 



(the set X^-'^ for clean dynamics) from the noisy data 
{xi}. Taking into account the above reasons it is useful to 
replace the global minimization problem ^ by N local 
minimization problems related to the nearest neighbor- 
hood X^^. The corresponding local functionals to be 
minimized are 

S 

\/n = l,...,N where {s : a;^ e X^^ or £ X„+i }(10) 

We can consider the minimization problem as a cer- 
tain approximation of (|SJ|. The global problem (|SJ) is 
equivalent to Eq. ^ with 2N unknown variables that 
should be found single-time. The problem IjlOl) is equiv- 
alent to a system of coupled equations that should be 
solved several times and as a result one gets an approxi- 
mate global solution. Writing Eq. (|10|) in the linear form 
i.e. calculating the zeros of corresponding derivatives and 
using Eq. one gets N linear equations as follows 

M„-<5X^f, = B„ Vn=l,...,7V, (11) 

where (^X^)^ = {Sxn, Sxk,Sxj,Sxn+i,Sxk+i,6xj+i,Xn)- 
The matrices M„ corresponding to (|10|l avoid the dis- 
advantages of i.e. they are not singular, their 
dimension is small and they do not substantially depend 
on the initial approximation of nearest neighbors. 
Matrix M„ for m = 1 is given by 
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Vector B„ has the form Bj = {0, 0, 0, 0, 0, 0, -G„(0)}. 
Note that this matrix in one-dimensional case is the same 
for LPNC method given in Ref. 1], but constraints and 
matrix M„ will essentially differ in higher dimensions for 
both methods. 



III. COMPARISON TO STANDARD LP 
METHODS 

Minimizations problems used in standard LP methods 
and in this method Eq. Q are not equivalent because 
in our case we do not have to estimate the Jacobi ma- 
trix at all. These differences in practice are as follows a) 
Eq. (PI is nonlinear against corrections Sxi. The approx- 



imation in this case means a corresponding linearization, 
b) For constraints in standard LP methods we do not 
know the exact values of Jacobi matrix A. The approx- 
imation means that Jacobi matrix A is estimated from 
noisy data. The LP methods look for subsequent correc- 
tions to noisy data by finding of a subspace tangent to 
an unknown attractor corresponding to the clean dynam- 
ics and projecting noisy data on this subspace. If one 
tries to estimate the position of the tangent subspace, 
what is equivalent to estimation of the Jacobi matrix A 
from noisy data, the range of the neighborhood should 
be larger than the magnitude of noise. Such a procedure 
should allow to distinguish between the dominating di- 
rection (connected with system dynamics) and random 
directions connected with a noise (see Figs.^ and^^)- 
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FIG. 1: The plot of the clean attractor (continuous line), 
noisy data connected to this attractor (black dots) as well as 
range of the nearest neighborhood taken under consideration 
to determine tangent subspace (rectangle). The level of noise 
for the case c) is so high that a linear approximation is not 
longer valid. 
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If the noise level is very high it is not possible to use the 
tangent subspace to find the attractor of the clean dy- 
namics since the range of the necessary nearest neighbor- 
hood would be very large and the linear approximation 
would be invalid (see Fig. ^) On the other hand, 

if we consider the minimalization problem we do not 
need to find the Jacobi matrix A but only to take into 
account the constraint equation (Q. Such an approach 
makes it possible to use a neighborhood smaller than the 
noise magnitude and in our approach the corresponding 
number of nearest neighbors equals to 2. Note that here 
we encounter a flow, so nearest neighbors searching is not 
so biased as local projection. To find two nearest neigh- 
bors to x„ we use the Delaunay triangulation |23 | 
and the method to find is given in Ref. jjy. Searching 
nearest neighbors by Delaunay triangulation is very time 
consuming. That is why we first look for the Nnn near- 
est neighbors by means of Euclidian distance minimiza- 
tion and then perform the Delaunay triangulation only 
on this nearest neighborhood. This approach is as fast as 
standard nearest neighbors searching. Accordingly to the 
needs, e.g. in online noise reduction of human voice, one 
can think about speeding up the method. As we men- 
tion in our approach we need only two nearest neighbors 
as opposed to standard LP methods, so looking only for 
the two nearest neighbors close in time would make the 
searching for closest neighbors very fast and the method 
robust against high non-stationarity. 

One can ask on the smallest sampling rate per cycle TZS 
which makes the method applicable. As we have checked 
the method works even for the rate TZS comparable with 
the m parameter used in Eq. Q but then the efficiency is 
smaller than for the standard LP method. The method 
works the best for more than 2 • m samples per cycle. We 
have taken the embedding window d ■ r used in nearest 
neighbor searching as long as one cycle. The method 
is robust against changing the number of taken nearest 
neighbors, as opposed to standard LP methods what will 
be seen in the next section. 



FIG. 2: The plot of the %R parameter on the gain parameter 



IV. NOISE REDUCTION: EXAMPLES AND 
APPLICATIONS 

Let us define the noise level Af 

Af=^^, (13) 

<^DATA 

where ct is a standard deviation of noise and a data is 
the standard deviation of data. The efficiency of noise 
reduction we calculate by means of the gain parameter 
which is defined as 

g = lOlogf^^\ (14) 

V ^red J 

where cr^ioise is the variance of added noise and ct^^^ 
is the variance of noise left after noise reduction. The 
gain parameter can be transformed into another pa- 
rameter: %R, which says how much noise is reduced: 
%R = (1 - <7f^Jal„,J ■ 100%. In Fig.Elthe dependence 
of the %R parameter on the basic parameter Q is pre- 
sented. The gain parameter Q is commonly used in the 
evaluation of the noise reduction because it gives more 
relevant information especially in the regime %R > 90%. 
We use the standard Lorenz model to evaluate the perfor- 
mance of noise reduction methods. The Lorenz system is 
described by a system of three coupled differential equa- 
tions 

y = px — y — xz (15) 
z — xy — (3z. 

The model is widely used for a description of Rayleigh- 
Benard instabilities in fluids ,23j and in quantum optics 
for laser dynamics We use standard parameters for 
this system, i.e. ^ = 10;p = 28;P = 8/3, for which the 
standard " butterfly" attractor can be observed. To verify 
our method in a real experiment we have performed anal- 
ysis of data generated by a nonlinear electronic circuit. 
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FIG. 3: The plot of the correlation dimension D2 versus the 
threshold e for the clean Lorenz system. 
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FIG. 4: The plot of the correlation dimension D2 versus the 
threshold e for the Lorenz system with noise M = 48%. 



The Chua circuit is one of the simplest electronic nonlin- 
ear system that exhibits chaotic behaviour psl . The 
nonlinearity comes from two paralel connected negative 
resistors which are realized by amplifiers with a corre- 
sponding feedback. The Chua circuit has been studied 
in the presence of noise added to the outcoming signal. 
The noise (white and Gaussian) has been generated by an 
electronic noise generator. The LPNCF scheme of noise 
reduction improves estimations of invariant parameters. 
Figs l3l5l present calculations of the correlation dimension 
D2 versus threshold e for the clean Lorenz system, the 
Lorenz system with noise and the latter parameter after 
noise reduction respectively. Using a standard procedure 
one looks for a plateau in an intermediate threshold re- 
gion. In fact one can observe that the plateau D2 ~ 2 
cannot be found at Fig. 0] corresponding to the noisy 
Lorenz system with the noise level TV = 48% but is well 
seen after the noise reduction with the LPNCF method 
(see Fig. El). 
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FIG. 5: The plot of the correlation dimension D2 versus 
the threshold e for the Lorenz system with noise after noise 
reduction. 



We have made a quantitative comparison between LP- 
NCF and GHKSS method. GHKSS method [H is the 
implementation of the standard LP approach that we 
think are optimal and most efficient. In case of both 
methods we use the same scheme of neighbors searching, 
i.e. the minimization of Eucleadian distance, in addi- 
tion for LPNCF method we apply then the Delaunay 
search which is not needed for GHKSS method. For 
large sampling rate we did always some time averaging 
with curvature correction that improve the gain in 
both methods. For all the calculations we use 8 iter- 
ations of both methods. The projection dimension for 
GHKSS 3- 12 and m for LPNCF method are in the range 
4 — 64. We can say that the complexity in programming 
of the both methods is comparable. The LPNCF method 
implemented at 2.5GHz computer is approximately 2-3 
times slower to be used for on-line noise reduction in 
voice and to receive a substantial improvement of the 
voice recognition. 

Figs. I6I8I shows the clean Chua attractor, with mea- 
surement noise TV = 46% and after noise reduction re- 
spectively. It is used the LPNCF approach for m = 3 — 12 
when sampling ratio was about 50 per full cycle. The ef- 
ficiency was g = 9.48 (%i? = 88.7%) when the GHKSS 
method did G = 9.33 {%R = 88.3%). 

We analyze the behavior of noise reduction by the LP- 
NCF and the GHKSS method against increasing sam- 
pling rate per cycle TZS (see Fig. It is clear that the 
efficiency of LPNCF method should increase for larger 
sampling rate TZS. In the figure one can see for large TZS 
that the LPNCF method is more efficient than GHKSS 
method while it is less efficient for small TZS. We compare 
the efficiency of these two methods for various noise levels 
TV (see Fig.lTUI). One can see that LPNCF method in this 
case is more efficient for high noise levels starting from 
30%. This is because for large noise levels it is difficult in 
GHKSS method to determine properly the tangent sub- 
space as it was explained in the previous section. The last 
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FIG. 6: The sampling corresponding to a clean trajectory in 
the Chua circuit (real experiment). 
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FIG. 7: The sampling received from the Chua circuit in 
the presence of a measurement noise M = 46%. Note the 
difference in scale. 
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FIG. 8: The sampling of Chua circuit received after the noise 
reduction applied to data presented at Fig. 13 
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FIG. 9; The efficiency of noise reduction by LPNCF and 
GHKSS method for different sampling rate TZS. Here the 
Lorenz system ^ was used {N = 48%, iV = 5000, N,,„ = 
20). 
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FIG. f 0: The efficiency of noise reduction by the LPNCF and 
GHKSS method for different noise level N ■ Here the Lorenz 
system |2j was used (7^5 = 66, iV = 5000, iV„„ = 20). 



comparison of these two methods shows the dependence 
on the number of regarded neighbors Nnn (see Fig. Illll . 
In the LPNCF method the parameter Nnn describes the 
number of neighbors that are used in preliminary search 
for candidates to the Delaunay procedure. It is shown 
in this figure that for small number of neighbors the LP- 
NCF method can be used without the loss of efficiency. 
Such a behavior is very useful for non-stationary data, 
when the large number of neighbors could not be found or 
when because of correlated noise one should omit neigh- 
bors close in time. The GHKSS method did some noise 
reduction for small Nnn because here most corrections 
come from the time averaging which alone made Q = 1 .2 
{%R = 80.9%). 

We applied successfully our method to noise reduction 
from human voice On Fig. E|we show a clean time 
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FIG. 11: The efficiency of noise reduction by the LPNCF 
and GHKSS method for various number of neighbors N„n. 
Here the Lorenz system jg^l was used (A/" = 48%, TZS = 66, 
N = 5000). 
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FIG. 12: The voice time series of sentence "Hello world, my 
name is Krzysztof Urbanowicz" . From the upper panel to the 
bottom there are clean signal then a series with decreasing 
measurement noise and a noisy signal after noise reduction 
respectively. 



series of the recorded sentence "Hello world, my name 
is Krzysztof Urbanowicz" (upper panel) , this time series 
with temporally decreasing measurement noise (middle 
panel) and after noise reduction (bottom panel). Noise 
reduction was made in windows of length N = 5000, with 
parameter m = 3 — 12. The voice was recorded with sam- 
pling 22050iJz what gives TZS 120. The embedding 
window d • T = 100, so it covers almost the whole cycle. 
Fig. US presents the efficiency of LPNCF and GHKSS 
methods. As it is suspected the LPNCF method did less 
for small noise levels (see for the comparison Fig. I10|l . 
Note that here we use larger number of neighbors than 
at Fig. ^1 i.e. Nnn — 60. For large noise level both 
methods work comparably. The gain of the noise reduc- 
tion for whole data set shown in Fig. E| is Q = 11 A 
i%R = 92.8%) for LPNCF and G = 11.7 {%R = 93.2%) 
for GHKSS. Such values of noise reduction improve sig- 
nificantly voice recognition for intermediate noise levels. 
After the performed noise reduction the background noise 
is not heard in the recorded signal. 
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FIG. 13: The efficiency of noise reduction of human voice by 
the LPNCF and GHKSS method for different noise level M 
{N = 5000). 
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In conclusion we developed the method of noise reduc- 
tion design for flows. It uses a nonlinear constraints that 
appear due to the continuous behavior of flows. To effi- 
ciently perform the noise reduction one needs to find only 
two nearest neighbors. The method is robust against 
input parameters estimation as well as for highly non- 
stationarity data. We applied with success the method 
for noise from human voice separating. 
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